Spatial Prediction with INLA

Introduction

This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation.

GMRF Background

Blangiardo and Cameletti (2015) section 6.1.

• Observations aggregated to disjoint areal regions indexed by \(i\).
• Each region has unique parameter \(\theta_{i}\).
• \(\mathcal{N}(i)\) is the set of indices of neighbors of region \(i\) and \(\mathcal{N}_{i} = |\mathcal{N}(i)|\) is the number of neighbor of region \(i\).
• Local Markov property: given \(\boldsymbol{\theta}_{\mathcal{N}(i)}\), \(\theta_{i}\) is independent of all other \(\theta_{j}\).

• Then the precision matrix \(\mathbf{Q}\) of \(\boldsymbol{\theta}\) is sparse because only neighbors have nonzero coprecisions.

• Besag-York-Molliè model:
  • Exchangeable random effects \(u_{i}\) with intrinsic conditional autoregressive (iCAR) structure.
  • \(u_{i}|\mathbf{u}_{-i} \sim \mathrm{N}\left(\mu_{i} + \sum_{j} a_{ij} (u_{j} - \mu_{j}) / \mathcal{N}_{i}, \sigma_{u}^{2} / \mathcal{N}_{i}\right)\)
  • (What are the \(a_{ij}?\)).
  • iCAR is an improper prior because covariance matrix not positive definite but this is ok for random effects.

SPDE Background

Find Lindgren et. al. (2011) and fill in motivation

• Spatial process \(\xi(\mathbf{s})\).

• Solve \((\kappa^{2} - \Delta)^{\alpha / 2}(\tau \xi(\mathbf{s})) = \mathcal{W}(\mathbf{s})\).
  • \(\kappa\) is a range parameter.
  • \(\Delta\) is the Laplacian.
  • \(\alpha\) is a smoothness parameter.
  • \(\tau\) a precision parameter.
  • Exact and sationary solution: \(\xi(\mathbf{s})\) is a Gaussian field with Matèrn covariance function.

• Finite element approximation $$.

Geostatistics Example

Toy dataset from Blangiardo and Cameletti (2015).

# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')

# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))

# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))

# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
toy_fit <- inla(
  y ~ -1 + intercept + f(spatial_field, model = spde),
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
toy_fit$summary.fixed

toy_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)

Bei Dataset and inlabru

Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005).

# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')

bei_corners <- vertices.owin(Window(bei))
bei_domain <- cbind(bei_corners$x, bei_corners$y)
bei_full_mesh <- inla.mesh.2d(cbind(bei$x, bei$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_full_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei, pch = '.', col = 'blue')

bei_full_spdf <- as.SpatialPoints.ppp(bei)
# CHECK PRIORS!
matern_full <- inla.spde2.pcmatern(bei_full_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = matern_full) + Intercept
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf)

lambda_full <- predict(bei_full_lgcp, pixels(bei_full_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_full)
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

plot(lambda_full['sd'])
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')

# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
n_quads <- 10
botleft <- cbind(runif(n_quads, 0, 950), runif(n_quads, 0, 450))
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
    cbind(
      botleft[r, 1] + c(0, 0, 50, 50),
      botleft[r, 2] + c(0, 50, 50, 0)
    )
  )})
bei_win <- do.call(
  union.owin,
  apply(botleft, 1, function(x){return(
    owin(x[1] + c(0, 50), x[2] + c(0, 50))
  )})
)
bei_hole <- bei[complement.owin(bei_win)]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)

plot(bei_hole, main = 'Observed Subregion', pch = '.', cols = 'black')

bei_hole_mesh <- inla.mesh.2d(cbind(bei_hole$x, bei_hole$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_hole_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
points(bei_hole, pch = '.', col = 'blue')

bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
# CHECK PRIORS!
matern_hole <- inla.spde2.pcmatern(bei_hole_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = matern_hole) + Intercept
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf)

lambda_hole <- predict(bei_hole_lgcp, pixels(bei_hole_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_hole)
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(lambda_hole['sd'])
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')

plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)

bei_samp_mesh <- inla.mesh.2d(cbind(bei_samp$x, bei_samp$y),
                              cutoff = 50, max.edge = c(50, 100),
                              loc.domain = bei_domain)
plot(bei_samp_mesh, asp = 1)
plot(Window(bei), border = 'blue', add = TRUE)
plot(Window(bei_samp), border = 'blue', add = TRUE)
points(bei_samp, pch = '.', col = 'blue')

bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
# CHECK PRIORS!
matern_samp <- inla.spde2.pcmatern(bei_samp_mesh,
                                   prior.sigma = c(0.1, 0.99),
                                   prior.range = c(5, 0.01))
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = matern_samp) + Intercept
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf)

lambda_samp <- predict(bei_samp_lgcp, pixels(bei_samp_mesh), ~ exp(mySmooth + Intercept))

# Plot posterior means and posterior sd.
plot(lambda_samp)
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

plot(lambda_samp['sd'])
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')

Bei Dataset with gridding

centers <- gridcenters(dilation(bei_window_full, 40), 40, 20)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
                     count = NA_integer_, area = NA_real_)

for(r in seq_len(nrow(bei_df))){
  bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
                         bei$x < bei_df$x[r] + dx &
                         bei$y >= bei_df$y[r] - dy &
                         bei$y < bei_df$y[r] + dy)
  bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(bei_df$count, nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')

# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(bei_full_mesh, as.matrix(bei_df[bei_df$area > 0, c('x', 'y')]))

# Initialize exponential covariance structure for SPDE.
full_spde <- inla.spde2.matern(mesh = bei_full_mesh, alpha = 2)

# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = full_spde$n.spde)
stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')

# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = bei_full_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))

# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')

# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)

# Fit the model with INLA.
bei_full_fit <- inla(
  count ~ -1 + intercept + f(spatial_field, model = full_spde),
  offset = larea, family = 'poisson',
  data = inla.stack.data(stacks),
  control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)

# Output posterior summaries.
bei_full_fit$summary.fixed

bei_full_fit$summary.hyperpar

# Extract posterior mean of latent spatial field.
index_pred <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- bei_full_fit$summary.linear.predictor[index_pred, 'mean']
post_sd <- bei_full_fit$summary.linear.predictor[index_pred, 'sd']

# Plot the posterior mean and SD of the latent spatial field.
plot(im(t(matrix(post_mean, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior Mean of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

plot(im(t(matrix(post_sd, nrow = length(unique(centers$x)), ncol = length(unique(centers$y)))), unique(centers$x), unique(centers$y)), main = 'Posterior SD of Spatial Field')
plot(bei_window_full, add = TRUE)
points(bei, pch = '.', col = 'black')

beihole_df <- data.frame(x = centers$x, y = centers$y,
                         count = NA_integer_, area = NA_real_)

for(r in seq_len(nrow(beihole_df))){
  beihole_df$count[r] <- sum(bei_hole$x >= beihole_df$x[r] - dx &
                             bei_hole$x < beihole_df$x[r] + dx &
                             bei_hole$y >= beihole_df$y[r] - dy &
                             bei_hole$y < beihole_df$y[r] + dy)
  beihole_df$area[r] <- area(Window(bei_hole)[owin(c(beihole_df$x[r] - dx, beihole_df$x[r] + dx), c(beihole_df$y[r] - dy, beihole_df$y[r] + dy))])
}

par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(beihole_df$count, nrow = length(unique(beihole_df$x)))), unique(beihole_df$x), unique(beihole_df$y), unitname = 'meters'), ncolcours = range(beihole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei, pch = '.', col = '#00000040')

Bei Dataset with Simpson et al. (2016) method

This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.

The pseudodata are constructed as follows.

• Let \(n\) be the size of the point pattern.
• Let \(p\) be the number of nodes of triangulation.
• Define the pseudodata as \(\mathbf{y} = (y_{1}, \dots, y_{p}, y_{p+1}, \dots, y_{p+n})'\) where \(y_{i} = 0\) for \(i = 1, \dots, p\) (the traingulation nodes) and \(y_{i} = 1\) for \(i = p+1, \dots, p+n\) (the observed points).
• Define \(\boldsymbol{\alpha} = (\alpha_{i}, \dots, \alpha_{p}, \alpha_{p+1}, \dots, \alpha_{p+n})'\) to encode the numerical integration scheme where \(\alpha_{i}\) is the numerical integration weight for \(i = 1, \dots, p\) (the traingulation nodes) and \(\alpha_{i} = 0\) for \(i = p+1, \dots, p+n\) (the observed points).

Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to the intensity at that their respective locations, observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to

\[ \prod_{i=1}^{p+n} \eta_{i}^{y_{i}} \exp(-\alpha_{i} \eta_{i}). \]

(Is there a missing \(\alpha_{i}\)?)

# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25

boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
mesh <- inla.mesh.create(
  boundary = boundary,
  refine = list(max.edge = MAX_EDGE_LENGTH)
)

plot(mesh, main = 'Mesh for Bei Site')
points(bei, pch = 19, cex = 0.25, col = 'red')

pts <- cbind(bei$x, bei$y)

# Contruct the SPDE A matrix for nodes and points.
nV <-mesh$n
nData <- dim(pts)[1]
LocationMatrix <- inla.mesh.project(mesh, pts)$A
IntegrationMatrix <- sparseMatrix(i = 1:nV, j = 1:nV, x = rep(1, nV))
ObservationMatrix <- rbind(IntegrationMatrix, LocationMatrix)

# Get the integration weights.
IntegrationWeights <- diag(inla.mesh.fem(mesh)$c0)
E_point_process <- c(IntegrationWeights, rep(0, nData))

# Create the psuedodata.
fake_data <- c(rep(0, nV), rep(1, nData))

# Fit model to full site.
spde <- inla.spde2.matern(mesh)
formula <- y ~ -1 + intercept + f(idx, model = spde) # No covariates.
data <- list(y = fake_data, idx = 1:(nV), intercept = rep(1, nV))
result_full <- inla(
  formula,
  data = data,
  family = 'poisson',
  control.predictor = list(A = ObservationMatrix),
  E = E_point_process,
  verbose = TRUE
)
summary(result_full)

Call:
c("inla(formula = formula, family = \"poisson\", data = data, E = E_point_process, ",  "    verbose = TRUE, control.predictor = list(A = ObservationMatrix))" )

Time used:
 Pre-processing    Running inla Post-processing           Total 
         0.2388         44.7314          0.0567         45.0269 

Fixed effects:
             mean     sd 0.025quant 0.5quant 0.975quant    mode   kld
intercept -5.8164 0.5817    -6.9936  -5.8135    -4.6551 -5.8083 4e-04

Random effects:
Name      Model
 idx   SPDE2 model 

Model hyperparameters:
                 mean     sd 0.025quant 0.5quant 0.975quant   mode
Theta1 for idx  2.597 0.0542      2.493    2.596      2.707  2.591
Theta2 for idx -4.243 0.1483     -4.554   -4.235     -3.973 -4.204

Expected number of effective parameters(std dev): 404.96(18.67)
Number of equivalent replicates : 14.20 

Marginal log-Likelihood:  -19395.22 

# Fit model to point site with holes.


# Fit model to quadrat-sampled site.

References

Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.

Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.

Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.

Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.